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Abstract 

We perform Random Phase Approximation (RPA) study of collective ex- 
citations in the bose-fermi mixed degenerate gas of Alkali-metal atoms at 
T = 0. The calculation is done by diagonalization in a model space composed 
of particle-hole type excitations from the ground state, the latter being ob- 
tained from the coupled Gross-Pitaevskii and Thomas- Fermi equations. We 
investigate strength distributions for different combinations of bose and fermi 
multipole (L) operators with L = 0,1,2,3. Transition densities and dynami- 
cal structure factors are calculated for collective excitations. Comparison with 
the sum rule prediction for the collective frequency is given. Time dependent 
behavior of the system after an external impulse is studied. 
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I. INTRODUCTION 



Collective excitation is one of the most prominent phenomena in quatum many-body 
systems such as liquid helium, electron gas, nucleus, etc. In the recently developed Bose- 
Einstein condensates of trapped atomic gases collective oscillations were the first of 

the dynamical phenomena discovered B . Collective oscillations are characterized by various 
quantum numbers related to, e.g., the shape of oscillation, or internal degrees of freedom of 
the constituents such as spin, isospin, etc. The oscillation frequency, the damping width, 
etc. depend on the interparticle interaction of the constituents, and thus provide a clue to 
unravel the dynamical correlation of the many-body system. 

Study of the properties of trapped neutral atoms has been extended to fermi systems Q 
where the occurrence of the Fermi degeneracy was observed, and also to the mixture of bose 
and fermi particles. The latter system with condensed bosons and degenerate fermions was 
recently realized experimentally [|5]||. This system is one of the typical example in which 
particles obeying different statistics are intermingled. Theoretical studies of the bose-fermi 
mixed system of cold atomic gases have been done for static properties [f^-|T2|l , for the phase 



diagram and phase separation |13|--[T5| , for the stability of the system |T6|J1y| and for collective 
excitations |JT~8| pO|l . 

In Ref. WE\ sum rule approach has been applied for collective excitations in the bose-fermi 
mixed system. Average excitation energies for the states with multipoles L = 0,1,2 were 
calculated for both in-phase and out-of-phase modes of the bose and the fermi particles, and 
the dependence on the bose-fermi interaction strength was studied. The sum rule approach 
is a powerful technique for collective states in quantum many-body systems, and has been 
successfully applied to the excitation of bose condensed systems, see Refs. p|j2l[1. It does 
not give, however, direct information on the eigenstates of the system, but rather an average 
behavior of the strength distribution for the adopted multipole excitation operator. For a 
more detailed investigation on the dynamical properties of the system, one would need a 
study of individual eigenstates. 
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In the present paper we extend the previous study |1| of collective excitations in the 
bose-fermi mixed system. We calculate full particle-hole type excitations by diagonalization 
of an RPA type matrix. The motive of this calculation is threefold: first of all, the calcu- 
lation allows us to study the excitation spectrum and the strength distribution in contrast 
to the sum rule method which focuses on the strength weighted average of the prescribed 
multipole operators. We study, for instance, the degree of collectivity of the excitations 
depending on the strength of boson-fermion interaction, and estimate the damping of the 
collective excitation albeit within the space of particle-hole excitations. Information on the 
wave function allows us to calculate observables such as the dynamical structure factor. The 
latter for bose-condensed systems is now becoming available experimentally by two-photon 



spectroscopy |22| and is being studied theoretically p3fl . Secondly, we compare the results 
with those obtained from the sum rule approach. This provides a check on the approxima- 
tion adopted in the calculation such as the model space truncation. The comparison also 
allows one to examine the structure of the low- and high-lying collective modes which was 
speculated in |Tj| through the mixing angle 9 of multipole operators. Third, we can predict 
the time-dependent behavior of the system for a given external perturbation. This process 
is actually the one that has been employed in the previous experimental study of collective 
excitations in BEC. An RPA study of the bose-fermi mixed system has recently been done 
in [pOfl , where the response for an external multipole field is formulated in the form of an 
integral equation. In the present calculation we approach the problem by diagonalizing the 
RPA matrix, by incorporating the discrete nature of the excitation in an isolated trapped 
system, and investigate, in particular, the properties of the strength distribution for a com- 
bination of bose and fermi operators and for various values of the bose-fermi interaction 
strengths. 

The content of the paper is as follows: in the next section we derive the RPA equation for 
the bose-fermi mixed system using the equation of motion for particle-hole type excitation 
operators. The single particle (hole) states are obtained in the mean-field calculation, i.e., 



by solving the coupled Gross-Pitaevskii and Thomas- Fermi equations. In Section [III 



we 



first briefly describe the parameters and the numerical procedure employed in the present 
calculation. We then turn to the detailed studies of the results obtained, including ground 
state density, single particle states, and strength distribution for each multipole. Comparison 
with the sum rule calculation is given. Transition densities and dynamical structure factors 
for some of the collective excitations are presented. We finally consider the time development 
of the system after an external multipole impulse on the system. Last section is devoted to 
summary and conclusion. 



II. FORMULATION 



We consider a dilute spin-polarized bose-fermi mixed system trapped in a spherically 
symmetric harmonic oscillator potential at T = 0. The system is described by the Hamilto- 
nian 



H = H + V b + V bf 



(1) 



with 



H = Jd 3 r^(r) 
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where and ip are the boson and the fermion field operators. We take, for simplicity, 
the mass m and the trapping frequency u to be the same for the boson and the fermion. 
The boson-boson and the boson-fermion interaction strengths for the pseudopotentials, g 
and h, are related to the s-wave scattering lengths a bb and a b f through g = Arrh 2 a bb /m, 
h = Airh 2 a b f /m, while the fermion-fermion interaction is omitted as we consider a polarized 
dilute system at low temperature. 

According to the Landau picture of quantum liquid, low lying excited states of the system 
may be described by quasiparticle excitations which have the structure of particle-hole (p- 
h) type excitations from the ground state. Small amplitude collective oscillations of the 



system, corresponding to the zero sound, are given by the coherent superposition of these 
p-h excitations and are well described by RPA for many fermion systems. The corresponding 
excitations in the bose condensed system have been formulated in the Bogoliubov-de Gennes 
type equation. The latter is essentially the same as the one obtained in RPA, and we consider 
below the two types of excitations on the same footing. 

In order to formulate RPA for the bose-fermi mixed system, we first determine the 
ground state in the mean field approximation. The obtained mean field provides us with 
single particle energies and wave functions. Let us expand the field operators in terms of a 
complete set of single particle wave functions as 



where bk and a a are the boson and the fermion annihilation operators in the single particle 
states specified by the wave functions (fik and ip a . They satisfy the standard commutation 
or anticommutation relations. The single particle states with quantum numbers k 01 a 
are determined by minimizing the energy expectation value. The trial wave function of 
the system is given by the product of N b bosons in a single state fc=o (r) and the Slater 
determinant of Nf different fermionic states ip a (r). From the stationary condition of energy 
for the variation of these wave functions with a given number of bosons JVj, and fermions 
Nf, we obtain the set of coupled equations, i.e., the Gross-Pitaevskii equation for the boson 
wave function 
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and the mean-field equation for the occupied fermion states 
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is the ground state density of fermions, //& is the boson chemical potential and e£ are the 
single particle energies for the fermion. Set of equations (f|) and (^) together with the 
fermion number condition 

"d 3 rp f (r) = N f (7) 

constitute a closed set of equations. The single particle wave functions are normalized to 
unity. In Eq. (f|) we approximate the factor ]V& — 1 by iV&. Once the single particle wave 
functions for the occupied states are obtained from Eqs. (§) and (^), the wave functions of 
the unoccupied states are calculated from a similar set of equations, i.e., 



k 2 1 

-— V 2 + -mujy + ^|0 o (f)| 2 + hp f (r) 



Mr) = 4Mr) (8) 
for bosons with k ^ 0, and the same equation (S) for fermions but for unoccupied states. 



Orthogonality of these wave functions to the occupied ones are automatically satisfied |24 

For a large number of particles with a smooth mean field potential, the Thomas-Fermi 
calculation provides a good approximation. Assuming this to be the case, we determine the 
fermionic ground state density p/ by using 

^-(6n 2 Pf (f)f 3 + l -rmoy + hN b \<p {f)\ 2 = e F (9) 

together with the equation (f|) instead of fully solving the coupled set of equations (U) and 
(H). The Fermi energy ep is determined by integrating the density Pf(r) so as to give the 
fermion number. Numerical consistency of this procedure will be shown in the next section. 
With pf and 0o so determined, the wave functions (fik (k ^ 0) and ij) a are obtained as 
described above. 

We assume below that the mean field is spherically symmetric, and the single-particle 
quantum numbers a and k involve standard radial and angular momentum quantum numbers 
(n,£,m). (Note that the spin quantum number is frozen.) We consider the case in which 
the fermionic ground state (Slater determinant) is m-closed and thus is consistent with the 
spherically symmetric mean field. We denote by p (h) the fermion single particle states 
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unoccupied (occupied) in the Slater determinant, i.e., those above (below) the Fermi energy 
ep. For the bosons at T = the occupied state is the lowest single-particle state k — 0. 

We now define the creation and annihilation multipole operators of the excited state 
\vLM) with angular momentum quantum numbers L, M and an additional quantum number 
v. 
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where X, Y, [/, V are the amplitudes to be determined in RPA, and (l p m p lh — \LM) is 
the Clebsh-Gordan coefficient. 

In RPA the operator Q is determined so as to satisfy the equation of motion 
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where the terms in the 1. h. s. which involve different combination of operators a), a, b\b 
from that in are neglected by the assumption of RPA. The state 



\uLM) = Qt LM \0) with Q vLM \0) = 



(14) 



is then an approximate eigenstate of the Hamiltonian together with the correlated ground 
state given by |0). The amplitudes satisfy the orthonormality condition 

E(^ L X p l L - Y£ h * L Y p i L ) + E(^l< - V£V&) = 5 W , (15) 

ph k^O 

Substituting Eq. ([U]) into Eq. (|i~3f) we obtain the eigenvalue equation in matrix form: 
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where submatrices A and B are given by 
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with 1 = In Eq. (0) R b nl (r) and R f nl (r) are the radial parts of the single particle 

wave functions and if), which are defined by n ; m (r) = R b nl {r)Yi m {6 , ip) and i/) n i m (f) = 
i?{ ; (r)Y"; m (^, (p) where Yi m (6,(p) are the spherical harmonics. We omitted C^iV^ -1 ) terms in 
this calculation. Actually we obtain the same set of equations if we replace by, bo in Eq. flT2|) 
with y/Nb, which is equivalent to the Bogoliubov-de Gennes equation for the bosonic system. 

The response of the system to an external field is represented by the transition matrix 
elements of the relevant operators which connect the ground state and an excited state. The 
external probe for a collective oscillation is assumed to be long-ranged and is given by the 
standard multipole operators, i.e., 

f _L r 2 ( L = 0) 

F L {r) = IV*H . (18) 

[ r L Y L0 (6) (L + 0) 

As the system is composed of two kinds of particles, we define operators 

F b = Jd 3 rft(r)F L (r)0(r), F[ = J d 3 r^(r)F L (r}^(r) (19) 
and their combination 

Ft = F[ ± Ft, (20) 



where F[ with r = f,b,+, — are resepectively called 'fermionic', 'bosonic', 'in-phase' and 
'out-of-phase' type operators. The transition amplitudes for these operators are calculated 
in RPA as: 

<0|F L >L0> = (0|[i^,gI M ]|0> = (0\[F[,Fl LO }\0) ± (0\[F b L , 2&JI0), (21) 
(0|[F/,F^ o ]|0) =^^^(W|L0)[(-l)^ L + (-l)^] Jdrr"R{ lh r x Rl lp , (22) 
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where A = 2 for L = and A = L for L ^ 0. The strength distribution for the multipole 
operator F[ is then given by 

£ - n„) | (0|F[|z/L0)| 2 . (r = /, 6, +, -) (24) 

V 

They measure the collectivity of the excited states with respect to the multipole operator 
which characterizes the shape (or volume for L = 0) oscillation of the system. 

More detailed information on the structure of the collective excitation may well by rep- 
resented by the transition densities 
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where p^(f),/S 6 (r) are the fermion and boson density operators: 

p/ (f) = ft(r)ij;(r), p b (r) = ft(r)(j)(r). (28) 



Recent development of the two-photon Bragg spectroscopy El allows us to obtain dynamical 
structure factors (response function) for the excitations in trapped atomic gases. They are 
related to the transition densities as 

si(q,n)^J2\^ L ( q )\ 2 8(n-n v ), (29) 

V 

KM = J drr 2 j L (qr/h) 5pl L (r), (30) 



where q denotes momentum transferred to the system and are spherical Bessel func- 

tions of order L. The long wavelength limit of the dynamical structure factor is proportional 
to the strength distribution for the multipole operators defined above. 



III. CALCULATION 



A. Numerical procedure 



We consider a boson-fermion mixture of potassium isotopes, i.e., 41 K (boson) and 40 K 
(fermion). We take the same value m = 0.649 x 10~ 25 kg for the boson and the fermion 
masses. The boson-boson interaction strength g is obtained from the scattering length 

— 15.13nm for 41 K- 41 K [26], while the boson-fermion interaction strength h is varied. (A 
negative value for the scattering length abf for 40 K- 41 K has been suggested in |26| , although 



not well established yet.) Both particles are assumed to be trapped in the spherical oscillator 



potential with the same trap frequency tu = 100Hz. The oscillator length £ = ^/h/muj is 
4.03/im for the adopted value of ujq. 

Number of bosons is fixed at Nb =1000, while that of fermions is calculated in the 
Thomas- Fermi approximation for a fixed Fermi energy ep = (6 x lMO) 1 / 3 /^. The latter is 
determined so that the last oscillator shell with a number of quanta iV = 2n + 1 is given by 
^Fermi = 17 at /i/g = 0. This gives the number of fermions Nf = 1140 at h/g = 0. Nf is 
dependent on the boson-fermion interaction strength, and it is assumed that all the subshells 
given by the quantum numbers (n, t) are m-closed, which gives, for instance, Nf = 1050 at 
h/g = 8 and Nf = 1227 at h/g = —6 for the given value of €f- 

The single-particle wave functions were obtained from the coupled Gross-Pitaevskii and 
Thomas-Fermi equations by expanding the wave functions in the harmonic oscillator basis 
for the given oscillator constant £. We included associated Laguerre functions up to 
n = 30 for fermions and n = 15 for bosons. 

The excitation energies and the wave functions for each multipole L are obtained by 
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diagonalization of the RPA matrix. To construct the particle-hole basis we included single 
particle states for bosons up to n = 7 for I = and n = 6 for I = 1,2,3, so that seven p-h 
configurations are included for each L — 0, 1, 2, 3. For fermions, four (particle) states above 
the Fermi level €f and up to four (hole) states below ep have been included for each I. The 
number of fermi p-h configurations is 240 for monopole, 464 for dipole, 672 for quadrupole 
and 1124 for octupole. The dimension of the RPA matrix is twice the sum of bose and fermi 
configurations for each multipole. The number of bose configurations is the same for the 
four multipoles considered in the present calculation. 



B. Static properties 

In the ground state the fermion density distribution is much broader than that for the 
boson due to Fermi pressure PJl7[|, although at large negative value of h/g the bose-fermi 



attraction tends to produce a larger overlap of the two kinds of particles as shown in Ref. |T7 
As the bose-fermi interaction becomes repulsive the fermions are squeezed out from the 
central part. This in turn causes a smaller overlap of bosons and fermions, and thus a 
relatively small net effect of the interaction on the binding energy. For the present choice 
of parameters the fermions are distributed outside of the boson 'core' around h/g ~ 7, 



forming a 'shell' like structure p|J^, p!3| , p^] . Note that this behavior would be changed if one 
adopts different values for Nf/Nb and g, which may be represented by a single parameter a 
introduced in Ref. 0. 

The condition for the validity of the Thomas-Fermi approximation used to obtain the 
above fermion density distribution may be expressed in terms of the local de Brogile wave 



length A(r) = h/p(r), where p(r) = y2m(ep — V e g) with fermion mean field potential 
Kff( r ) = \vnijj\r 1 + hp h {r). With this quantity, the condition becomes [f2~7j 

d\{r) 



f(r) 



< 1. (31) 



dr 

At h/g = 8 where the fermionic potential may have a most pronounced structure, the value 
of f(r) is ~ 10 -2 except around turning points. At the turning points, however, the fermion 
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density produces only a negligible influence on the mean field potential. The validity of the 
present Thomas-Fermi calculation may also be checked by comparing the fermion density 



h/g — —6 ~ 8 shows that they agree within the order of 10~ 2 over the entire radial range, 
suggesting the consistency of the present calculation. 

In Fig. [I] we show the fermion single-particle energies measured from Fermi energy, 
e» — £Fi against orbital angular momentum £ for three values of the interaction parameter, 
h/g = 5.0, 1.0, —3.0. The Fermi energy €f is denoted by the horizontal line at zero energy. 
Note that e a are given by those of the harmonic oscillator at h/g = 0, e n ^ = (2n+£+3/2)huo, 
as we have no direct interaction among fermions. One should note also that the yrast state, 
i.e., the lowest state for each angular momentum, has no radial node, and the yrare one, the 
second lowest state, has one radial node, etc. The nodal structure of the particle and hole 
states around the Fermi surface is responsible for the multipole strength distribution of low 
energy particle-hole excitations. The figure shows that the states with low orbital angular 
momentum are much influenced by the bose-fermi interaction, while those with high angular 
momentum are almost insensitive to the values of h/g. This is because the additional fermion 
potential due to the bose-fermi interaction vanishes outside the boson density distribution. 
For instance, since the yrast single particle wave functions with angular momentum £ are 
peaked around v££, those fermion states with £ > (_Rb/£) 2 , Rb being a typical edge radius 
of the boson distribution, would not be much influenced by the bose-fermi interaction. (In 
the present case Rb — 3£, and the above relation gives £ > 9.) 

C. Energy weighted moments and comparison with the sum rule calculation 

RPA calculation provides approximate eigenvalues and wave functions for all the indi- 
vidual eigenstates and is useful to study details of the dynamics of the system. If one is 
interested in the gross behavior of the strength distribution or an approximate frequency of 
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the oscillation, one may rather consider the p-th energy weighted moments of the strength 
distribution, 



m 



,(L,T)=J2(hn„T\(0\F»\ 2 , (32) 



where hVL u is the excitation energy of the state \v). These moments can be expressed as 
ground-state expectation values of multiple commutators of F[ with the Hamiltonian |29| , |30| . 
It is known that this relation, sum rule, is conserved in RPA for some of the moments, i.e., the 
sum in the r.h.s. of Eq. (|3~2D obtained from RPA eigenstates coincides with the expectation 



value of the commutator in the HF ground state p8|)P9|. Thus a comparison of the two 
would provide a criterion on the consistency of the RPA calculation. In particular, the first 
moment (energy weighted sum, EWS) rri\ is calculated from the double commutator of F£ 
and the Hamiltonian as 

mi {L,T) = \{[Fl,[H,Fl]]) Q , (33) 



where ( )o denotes a ground state expectation value p§| . We have checked that the sum 
rule is satisfied within 1% for all the multipole modes. 

One can estimate the average frequency of the collective oscillation based on sum rules. 
There are several ways to define the average frequency depending on which part of the 
strength distribution is emphasized. In Ref. [18| the ratio of the third and the first moments 




(u = J'-!* (34) 



was studied based on sum rules. This definition is advantageous as it can be used to test 
the validity of the RPA calculation as mentioned earlier. In later discussions we consider 
also the average frequency calculated by 

u av = — (35) 
m 

which has more weight in the low-frequency strength compared with tD. The two frequencies 
should coincide when a single collective state exhausts the strength. We discuss later also 
the width of the strength distribution defined by 
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(36) 



In Fig. H we show the average frequency Eq. ( p4|) for L = 0,1,2 and r = +, — given 
by the sum rule ([33|) and the one directly calculated in RPA. We find a good agreement of 
the two calculations in a wide range of the bose-fermi interaction parameter h for a fixed g. 
There is a discrepancy at very large values of \h\, especially in the strongly attractive case. 
In the latter case, an induced instability of the ground state due to the bose-fermi attraction 



may be close |]I6 |, which suggests that a larger configuration space may be required to satisfy 
the sum rule. In fact, by increasing the number of particle-hole states, we obtained a better 
agreement. 



D. Distribution of multipole strengths 

Now we show the distributions of the multipole strengths. Figures |3],|7|, |lT||Tl"| respectively 
show the strength distributions for L =0,1,2 and 3, for either the fermionic/bosonic or the 
in-phase/out-of-phase operators. 



1. monopole 

In Fig. ||we show the monopole strength distribution, |(0|Fq |z/)| 2 against Q u , for r = f,b 
and for three values of h/g. Here one expects a volume oscillation of boson and/or fermion 
densities. For h = and for a large number of particles, the bosonic monopole oscillation 
will be located around \fhuJo as given by the collisionless hydrodynamics 0,^TJ, while the 
fermionic one is concentrated at 2ujq. The latter property is due to the almost degenerate 
values of the relevant particle-hole energies contributing to the monopole oscillation. The 
present calculation show that this situation persists even for large values of h/g, although 
the bosonic frequency is slightly shifted. The fermionic strength is distributed over a few 
states around 2u , showing that the induced fermi-fermi interaction via bosons is not strong 
enough to make a coherent superposition of the fermion particle-hole states. The total 
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monopole strength for fermion is by an order of magnitude larger than that for boson, even 
though the number of fermions involved in the excitation is smaller because of the Pauli 
principle, e.g. mi = 2590fru;o£ 4 for fermion and mi = 419hu>o^ 4 for boson at h/g — 1. This 
is because of the broader density distribution for fermions. 

The strength distribution may suggest that the fermions and bosons are moving indepen- 
dently even for rather large values of h/g. This is not necessarily the case, however, as one 
may see in Fig. ||, where we plot the dynamical structure factor for the three cases, i.e., the 
states at Q/ujq = 2.03 for h/g = 5.0, at 2.01 for h/g = 1.0 and at 1.98 for h/g =-3.0, which 
carry the largest strengths for the in-phase monopole operator. Although the fermionic 
strength is far larger than the bosonic one for this state, the mixture of the bosonic compo- 
nent is not small and is peaked at larger values of the momentum transfer q as shown by 
the dashed lines. The latter shows that the bosons oscillate in the inner region although not 
quite recognizable as far as one studies only the F distribution. 



In Refs. [|16|,|T8j we suggested that an instablity towards collapse may occur at large 
negative values of h/g and may be signaled by the lowering of uo. In the present calculation 
with smaller number of particles this is not apparent in Fig. |2|. In Fig. [5] we show the strength 
distribution at h/g = —6.65. This value is close to the critical value of instability around 
—6.7 estimated from the Maimer's condition [f7,16|| and around —6.8 from the condition by 



Roth and Feldmeier [ljj. We find a lowering of a single state which carries 13.9% of the 
EWS. The average frequency for this case is uo = 1.92co>o appreciably lower than the value 
~ 2uq for h/g = —5 ~ 5. At more negative values of h/g we could not find a stable ground 
state. 

Let us study the character of the low-lying excited states by considering the response to 
the probe 

F(9) = F + cos6 + F~ sinO (37) 

parametrized by e. In Ref. the angle e was determined by minimizing the average 
frequency u for F(e) so as to find the character of the probe which favors the low-lying 
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states. Once the value of the parameter, ^min is determined for a given h/g, one may consider 
an operator 'perpendicular' to F(9 min ), i.e., F± = F + sin 9 min — F~ cos 9 min , which may favor 
the high-lying state. Alternatively, one may maximize the EWS within the model space 
to determine the value # max of the operator to characterize the high-lying states. Figure |6| 
shows the strength distribution for these three types of operators at h/g = 5.0. The value 
#min = — 0.127T suggests that, at this highly repulsive value of the interaction, a slightly in- 
phase type oscillation is favored in order to avoid an overlap of bosons and fermions, see the 
density distribution in Ref f^]8], |i3l , |T7ll . We note that the distribution for the 'perpendicular' 



operator F± is concentrated and is similar to the one for F(9 max ). The latter, however, 
is an almost bosonic type operator and suggest that the strength distribution alone is not 
sufficient to characterize the structure of high-lying states. A similar study has been made 
for a very attractive case h/g = —6.65. In this case we obtain # min = — 0.017T, i.e., almost 
in-phase, and the strength distribution is similar to Fig. |5|. Thus the in-phase character of 
the low-lying mode which favors the overlap of the two kinds of particles is clearly seen in 
the determination of F(9 min ). 

2. dipole 

Dipole strength distributions for the in-phase and the out-of-phase operators are shown 
in Fig. |7|. The strong peak in the in-phase strength distribution corresponds to the center- 
of-mass oscillation of the whole system. For many-particle systems confined in a common 
oscillator potential the center-of-mass motion is decoupled from other (intrinsic) degrees 
of freedom of the system, resulting in the oscillation with the same frequency u . This 
relation can be represented by the commutation relation of the dipole operator and the 



Hamiltonian, and should hold also within RPA |[28|| . In the present numerical calculation, 
however, the state at Q = does not exhaust the whole strength for the in-phase oscillation; 
the largest deviation occurs at h/g — — 1 having 80% of the EWS. This is because the single- 
particle model space in our numerical calculation is not sufficient to completely decouple 
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the center of mass motion, especially at strong coupling cases where the deviation of the 
single particle potential from the oscillator becomes large. The mi sum rule is, nevertheless, 
almost satisfied as mentioned earlier, and is determined by the total number of bose and 
fermi particles, where the latter number is dependent on the value of h/g. Smaller value 
of the sum rule percentage at h/g = — 1 is mainly due to the appearance of the almost 
degenerate state (Au/uj ~ 1CT 5 ) at this energy. 

All other dipole oscillations should be orthogonal to the center-of-mass motion and thus 
are out-of-phase type in character. One may note that at repulsive values of h the dominant 
part of the out-of-phase strength is located below ujq, although the strength carried by each 
individual state is not very large. As discussed in Ref . |1| , this may reflect the ground state 



density distribution which favors the out-of-phase oscillation by making the overlap of boson 
and fermion distributions smaller. Figure § shows transition densities for the out-of-phase 
type oscillations at h/g = 5.0,1.0,-3.0 which carry the largest strengths; Q/u = 0.700 
for h/g = 5.0, 0.956 for h/g = 1.0 and 1.14 for h/g =-3.0. Bosons and fermions move in 
the opposite directions so as to make the overlap smaller around the surface of the boson 
distribution. One can observe that the bosonic transition density is large and robust, while 
the fermionic one is feeble and is spread over the whole system. This situation may be 
analogous to the soft dipole mode speculated in neutron-rich nuclei wherein the protons 
oscillate almost free in the sea of neutrons f3lf . Dynamical structure factors corresponding 
to these states are shown in Fig. ^. Here the bosonic and the fermionic responses are 
shown together with the one for a hypothetical out-of-phase type probe. By changing the 
momentum transfer one would find a structure corresponding to the oscillation of the two 
kinds of particles. 

It was suggested in ]7|,[19| that for a very strong repulsion between bosons and fermions 
with sufficient number of particles, the system may become unstable towards a phase separa- 
tion of the two kinds of particles. Out-of-phase dipole strength distribution in Fig. [7] indeed 
show a softening of the strength distribution at large value of h/g. We cannot conclude from 
the present calculation, however, if this tendency is related to the mentioned instability. 
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3. quadrupole 



One finds from the quadrupole strength distribution in Fig. 10 that the fermionic strength 
is split into low- and high-energy parts, while the bosonic one is concentrated. The higher 
strengths for the fermions are due to the 2huo excitation of the particle from the occupied 
single-particle states and are similar in character to the bosonic excitation. In contrast, 
the lower fermion strengths come from the matrix elements of the quadrupole operator F2 
which re-orient the single particle states within the same oscillator shell of N = In + £ 
around the Fermi surface. This transition normally involves a change of nodes by one, see 
Fig. |l], and the corresponding strengths are smaller than the higher ones. One may note 
that the quadrupole bosonic strength is located slightly higher than v^2^o as expected for 
the collective oscillation in the large N limit || , which may be due to a rather small number 
of particles in the present calculation. 

We note that the strength distribution tends to be fragmented at large (repulsive) value 
of the interaction parameter h/g as seen in Fig. |H]. This is because the average potential 
deviates appreciably from the harmonic oscillator potential, giving rise to a dispersion in 
the fermion particle-hole energies, see Fig. [I]. Strong fermion-boson interaction would then 
be responsible to scatter the bosonic strengths, too. 



4. octupole 



Octupole strength distribution is shown in Fig. [11]. Here again the fermionic strength 
is split into hu and 3hu>o regions due to the character of the multipole operator F 3 . We 
note also that the frequency of the bosonic oscillation is much larger than the hydrodynamic 
value v^u'o pl| . 

We find a striking change in the fermionic strength distributions depending on the sign 
of the bose-fermi interaction. While for a strongly repulsive case {h/g = 5.0) the region 
between the low- and high-lying states is almost filled up by small strengths, the one for a 
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strongly attractive case (h/g = —3.0) shows a large gap in the strength distribution. This 
may be traced back to the structure of the single-particle states in Fig. |I[ The large strengths 
around 3u>o are due to high orbital angular momentum states with £ p — = 3 which stay 
robust against h/g as discussed earlier. Smaller strengths due to low £ p and £h states are, 
on the other hand, sensitive to the interaction, and the correponding particle-hole energies 
are below (for h/g > 0) or above (for h/g < 0) the unperturbed value 3ujq. This may be 
understood by noting that the particles passing through the center (small I orbit) feel the 
shallow (for h/g > 0) or deep (for h/g < 0) potentials and are thus more easy or hard to 
excite. For the strengths around ujq different combinations of orbits are involved, and the 
strengths are broadened by the deviation from the oscillator potential. 

E. Time evolution 

To study the collective oscillation of the cold atomic gases, the current experiments on 
BEC exert a time-dependent field on the system and then observe the time development of 
the shape and the size of the condensate 0. This procedure generally excites a number of 
normal modes, and one may in principle be able to resolve each mode by performing Fourier 
transform as far as the time duration is long enough. To simulate the situation we consider 
a time evolution of the system after one applies an external weak pulse of the step-function 
type, i.e., AV oc F[ for —At < t < 0. The calculation is performed within the lowest order 
perturbation theory. (At is chosen as At = lCT 2 ^ 1 .) 

1. monopole 

First we consider a perturbation of monopole type 

AV = muo Auj VAitFq (38) 

with Au <C ujq which is applied to the ground state of the system for a short period At. The 
effective frequency of the external oscillator potential for fermions and/or bosons is then 
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changed into ujq ± Aoo depending on the choice of r. Time dependence of the rms radius 



of bosons and fermions is plotted in Fig. |P2 at h/g — —3 for the fermionic external field, 



t = f. Frequency of the short period oscillation is related to the average frequency of the 
strength distribution, e.g., u av = 1.99u;o- That of the long period one, on the other hand, 
reflects the width of the distribution. For instance, the slow decrease of the envelope for 



(rj) would imply a very narrow width of the strength distribution. The bottom of Fig. |12| 
shows the Fourier transform of (r^)(t) — (^/)o which recovers the sharp peak structure of the 
strength distribution in Fig ^. Time evolution of the bosonic radius shows a rather regular 
modulation, and its beat frequency is estimated about 0.125co>o, which is consistent with the 
half width a/2 = 0.127w . 

2. dipole 

For the dipole case we take the perturbing potential 

Ay = - mWo 2 A^-F[, (39) 

where Az(-C £) measures the shift of the centers of the oscillator potentials felt by fermions 
and/or bosons. In Fig. [H| we show the time dependence of the center-of-masses of bosons and 
fermions for the out-of-phase external dipole impulse, r = — , at h/g = 5. The frequency 
of the oscillation is consistent with the average frequency 0.90cg>o for the dipole strength 
distribution. The amplitude of the oscillation shows an irregular behavior, suggesting that 
the strength distribution does not have a simple structure. The lowest part of the figure 
shows the Fourier transform of the relative distance of the fermion and boson center-of-mass 
positions, Nf(zf) —Nb{zb), which is sufficient to recover the behavior of the original strength 
distribution given in Fig. [7|. 
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IV. SUMMARY AND CONCLUSION 



In the present paper we performed an RPA calculation of the polarized bose-fermi mixed 
system of alkali-metal gases at zero temperature. We solved a coupled Gross-Pitaevskii- 
Thomas-Fermi equations to obtain the ground state for several values of the boson-fermion 
interaction strength with a fixed value of the boson-boson interaction. Single particle states 
of bose and fermi particles are calculated based on the mean field produced by the obtained 
ground state density. The density distribution constructed from the occupied single-particle 
orbits agrees well with the original one, suggesting the self-consistency of the calculation. We 
calculated and diagonalized the RPA matrix to obtain excitation energies and wave functions 
for multipoles L — 0, 1, 2, 3. We then calculated strength distributions for bosonic/fermionic 
and in-phase/out-of-phase multipole operators, and also transition densities and dynamical 
structure factors for some of the collective states. 

We first calculated energy weighted moments of the strengths to study average behavior 
of the strength distribution. Comparison of the first moments with the sum rule predic- 
tions for L — 0, 1, 2 shows that the RPA configuration space is sufficient to within 1%. A 
glance at the monopole distribution suggests that fermions and bosons are moving rather 
independently for the adopted parameters , which turns out not to be the case if we study 
the dynamical structure factors. We find, in fact, that the long wavelength oscillation of 
fermions in the surface is coupled to the internal short-wavelength oscillation of bosons. 
For a strong attractive boson-fermion interaction, the calculation suggests a softening of in- 
phase monopole oscillation. We studied also the structure of the low-lying out-of-phase type 
dipole modes. Lowering of the energy is related to the ground state density distribution of 
bosons and fermions which favors the out-of-phase oscillation. Transition densities for these 
modes suggest that the bosons oscillate on their way through the cloud of fermions, which 
is analogous to the soft dipole mode discussed in neutron-rich nuclei. We also calculated the 
strength distributions for quadrupole and octupole modes, and studied, in particular, the 
origin of the fragmentation of fermionic strengths. Finally we considered time dependent 
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behavior of the trapped bose-fermi system after an impulse of a multipole external field. 
By Fourier transforming the time- dependent oscillating behavior of the system, we could 
recover the gross structure of the strength distribution. 
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FIGURES 

FIG. 1. Fermion single particle energies for h/g = 5.0,1.0,-3.0 against orbital angular mo- 
mentum quantum number. Energies are measured from the Fermi energy. Eigenstates which are 
degenerate at h/g = are connected by solid lines. 

FIG. 2. Average frequencies of the collective excitations in the sum rule formalism. Average 
frequencies u, Eq. (|34|), in units of ujq is plotted against interaction strength ratio h/g with fixed 
g. Solid and dotted lines respectively show the in- and out-of-phase oscillations. For a comparison, 
average frequencies calculated from the ground state expectation values of double commutators , 
see Ref. are plotted by dashed (in-phase) and dot-dashed (out-of-phase) lines. 

FIG. 3. Monopole strength distribution for fermionic and bosonic operators, Fq and Fq, at 
h/g =5.0, 1.0 and -3.0. Strengths in units of £ 4 are plotted against excitation energy. 

FIG. 4. Dynamical structure factors for the states in which in-phase monopole strengths are 
the largest. These states correspond to the ones at Q/ujq = 2.03 for h/g = 5.0, at 2.01 for h/g = 1.0 
and at 1.98 for h/g =-3.0. Dotted and dashed lines correspond to the fermion and boson transition 
densities, 5p{,Q , while solid lines are for the in-phase one, 5p^ = Spl + 5p b l/0 . The abscissa is q£/h, 
where q denotes the momentum transferred to the system by external probes. 

FIG. 5. Strength distribution for the in-phase monopole operator Fq at h/g =-6.65 plotted 
against U/ujq. 

FIG. 6. Strength distribution for the ^-dependent monopole operator Fo(9) at h/g = 5.0. 
Upper figure shows the strength distribution for the operator F{6 ia \ n = — 0.127r) which was de- 
termined so as to minimize the average frequencey u. Middle figure is the one for the operator 
F± = F{6 = 0.387r) which is perpendicular to the above operator. Lower figure is the strength 
distribution for the operator F{6 rafai = 0.277r) which was determined to maximize the average 
energy. 
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FIG. 7. Dipole strength distribution for in- and out-of-phase operators, and F 1 . Strengths 
are measured in units of £ 2 . 

FIG. 8. Transition densities for the dipole oscillations with largest strengths of the out-of-phase 
type. Excitation energies of these states are Q,/lvq =0.700 for h/g =5.0, 0.956 for h/g =1.0 and 
1.14 for h/g =-3.0. Solid and dotted lines respectively show the fermionic and bosonic transition 

f b ^ 

densities, 5p u \ , in units of £ . 

FIG. 9. Dynamical structure factors for the same states as given in Fig. [8|. Dotted and dashed 
lines are those for the fermionic and bosonic transition densities while solid lines are those for the 
out-of-phase one, 5p~ x = 5p^ — 5p b yl . 

FIG. 10. Quadrupole strength distributions for fermionic and bosonic operators, i 7 / and 
at h/g =5.0, 1.0 and -3.0. 

FIG. 11. Octupole strength distributions for fermionic and bosonic operators, Fg and F%, at 
h/g =5.0, 1.0 and -3.0. 

FIG. 12. Time dependent oscillation of the bose-fermi system at h/g =-3.0 after an external 
impulse of fermionic monopole-type, Eq. p^). Time dependences of the root-mean-square radii 
(in units of £) of fermions (upper figure) and of bosons (middle figure) are shown. Parameters of 
the impulse are: Aw/wo = 0.01 and uoAt = 0.01. Differential equation in time is solved with the 
time- mesh of 0. 05^(7 1 up to t max = 100/uq. The lowest figure shows the Fourier transform (in an 
arbitrary unit) of the fermionic rms radius deviation from the ground state, (r 2 j){t) — (r'j)o, which 
corresponds to the fermionic strength distribution given in Fig. || at h/g = —3.0. 
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FIG. 13. Time dependence of the oscillations of fermions and bosons at h/g = 5.0 for an 
out-of-phase type external dipole impulse, Eq. (j39|). Upper and middle figures show the fermion 
and boson center-of-mass in the unit of £ against elapsed time in the unit of uJq 1 . Parameters are 
given by Az/£ = 0.01 and ajgAt = 0.01. Calculation is performed with time-mesh of 0.05UJQ 1 and 
up to t max = IOO/cjo- The lowest figure shows the Fourier transform (in an arbitrary unit) of the 
relative distance of fermions and bosons, Nf(zf) — Nb(zb), with an energy resolution of ~ 10~ 2 . 
The figure corresponds to the strength distribution of the out-of-phase type in Fig. [j] at h/g = 5.0. 



28 



